// STL
#include <cmath>                               // std::sqrt()
#include <vector>                              // std::vector<>

// Boost
#include <boost/math/distributions/normal.hpp> // boost::math::normal, ::cdf(), ::quantile()

// stats++
#include "statsxx/postprocess/ROC.hpp"         // ROC_pt


//
// DESC: Calculate (explicitly) the binormal ROC curve.
//
// INPUT:
//
//     double                             a   : binormal coefficient
//     double                             b   : binormal coefficient
//
//     unsigned int                       n   : number of points to calculate binormal ROC curve over (default of 100).
//                                              NOTE: The actual number of points stored in the ROC curve will actually be two less, because the end points are discarded.
//
// OUTPUT:
//
//     std::vector<ROC_pt>                ROC : ROC curve
//
// -----
//
// NOTE: Scores are not calculated for the (returned) ROC curve.
//
// NOTE: In an earlier version, ROC scores were first mapped directly to updated (FPR,TPR) points obtained by Deming regression, and then spline interpolation was done on the line in Phiinv space.
//
// NOTE: ... Now, map_scores() is used as a standard for all fittings.
//
// -----
//
// NOTE: See notes in fit_binormal() for relevant formulas and explanations.
//
inline std::vector<ROC_pt> binormal_to_ROC(
                                           const double        a,
                                           const double        b,
                                           // -----
                                           const unsigned int  n
                                           )
{
    std::vector<ROC_pt> ROC;

    //=========================================================
    // INITIALIZE
    //=========================================================

    boost::math::normal dist;

    //=========================================================
    // CALCULATE BINORMAL CURVE
    //=========================================================

    // NOTE: See the NOTEs in ROC.hpp, related to / justifying plotting uniform in FPR.

    //---------------------------------------------------------
    // CALCULATE FPR AND TPR AND Phiinv_FPR AND Phiinv_TPR
    //---------------------------------------------------------
    std::vector<double> FPR;
    std::vector<double> TPR;
    std::vector<double> Phiinv_FPR;
    std::vector<double> Phiinv_TPR;

    for( unsigned int i = 1; i < (n-1); ++i )
    {
        double _FPR        = static_cast<double>(i)/static_cast<double>(n-1);
        double _Phiinv_FPR = boost::math::quantile( dist, _FPR );

        double _Phiinv_TPR = a + b*_Phiinv_FPR;
        double _TPR = boost::math::cdf( dist, _Phiinv_TPR );

        FPR.push_back(_FPR);
        TPR.push_back(_TPR);
        Phiinv_FPR.push_back(_Phiinv_FPR);
        Phiinv_TPR.push_back(_Phiinv_TPR);
    }

    //=========================================================
    // CONSTRUCT ROC CURVE
    //=========================================================

    for( auto i = 0; i < FPR.size(); ++i )
    {
        ROC_pt pt;
        pt.FPR = FPR[i];
        pt.TPR = TPR[i];

        ROC.push_back(pt);
    }

    // ADD IN THE ASYMPTOTIC POINTS
    ROC_pt pt;

    pt.FPR = 0.;
    pt.TPR = 0.;

    ROC.insert(ROC.begin(), pt);

    pt.FPR = 1.;
    pt.TPR = 1.;

    ROC.push_back(pt);

    //=========================================================
    // RETURN
    //=========================================================

    return ROC;
}
